In [1]:
import DSP
using PyPlot
In [2]:
include("../juwvid.jl")
Out[2]:
In [3]:
### Using LIGO data, provided by https://losc.ligo.org/events/GW150914/
data=readdlm("fig1-observed-H.txt"); # Hanford: Get from https://ligo.caltech.edu/WA
t=data[:,1]
y=data[:,2];
data=readdlm("fig1-observed-L.txt"); # Livingston: Get from https://ligo.caltech.edu/LA
t2=data[:,1]
y2=data[:,2];
In [4]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1,aspect=0.05)
PyPlot.plot(t,y)
PyPlot.xlabel("time [s]")
PyPlot.xlim(0.25,0.46)
PyPlot.ylim(-1.5,1.5)
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2,aspect=0.05)
PyPlot.plot(t,y2)
PyPlot.xlim(0.25,0.46)
PyPlot.ylim(-1.5,1.5)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
Out[4]:
In [5]:
# compute frequency indices corresponding to 0-500 Hz
dt=t[2]-t[1]
nufft=1024
juwutils.frequency_to_index([0.0,500.0], dt, length(y),nufft)
Out[5]:
In [6]:
# define the frequency range
fin=collect(linspace(1.0,211.0,nufft));
# compute the analytic signals
z=DSP.Util.hilbert(y);
z2=DSP.Util.hilbert(y2);
In [7]:
tfrst=stft.tfrstft(y);
tfrst2=stft.tfrstft(y2);
In [8]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(abs.(tfrst).*abs.(tfrst),t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
PyPlot.ylim(0,500)
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(abs.(tfrst2).*abs.(tfrst2),t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
PyPlot.ylim(0,500)
Out[8]:
In [9]:
include("../cohenclass.jl")
tfrwv=cohenclass.tfrwv(z,NaN,NaN,fin,NaN,0,"nufft");
tfrwv2=cohenclass.tfrwv(z2,NaN,NaN,fin,NaN,0,"nufft");
In [10]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(tfrwv,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(tfrwv2,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
Out[10]:
In [11]:
tfrpwv=cohenclass.tfrpwv(z,NaN,NaN,fin,NaN,NaN,0,"nufft");
tfrpwv2=cohenclass.tfrpwv(z2,NaN,NaN,fin,NaN,NaN,0,"nufft");
In [12]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(tfrpwv,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(tfrpwv2,dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
Out[12]:
In [18]:
#alias free sm
afwv=smethod.tfrsm(y,NaN,4,NaN,2)
afwv2=smethod.tfrsm(y2,NaN,4,NaN,2)
#L wigner distribution (L=2)
tfrlw=lwigner.tfrlw2L(afwv,4);
tfrlw2=lwigner.tfrlw2L(afwv2,4);
In [19]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(tfrlw,t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
PyPlot.ylim(0,500)
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(tfrlw2,t[2]-t[1],t[1],t[end],NaN,NaN,11.0)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
PyPlot.ylim(0,500)
Out[19]:
In [20]:
tfrpo=polywv.tfrpowv(y,NaN,NaN,fin,8,8);
tfrpo2=polywv.tfrpowv(y2,NaN,NaN,fin,8,8);
In [21]:
fig=PyPlot.figure(figsize=(10,5))
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(log10.(complex(tfrpo)),dt,t[1],t[end],fin[1],fin[end],0.7,"CMRmap",4,10)
PyPlot.colorbar(a,shrink=0.3)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(log10.(complex(tfrpo2)),dt,t[1],t[end],fin[1],fin[end],0.7,"CMRmap",4,10)
PyPlot.colorbar(a,shrink=0.3)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
Out[21]:
In [22]:
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.wtfrshow(tfrpo./maximum(tfrpo,1),dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.xlabel("time [s]")
PyPlot.ylabel("frequency [Hz]")
PyPlot.title("Hanford")
ax = fig[:add_subplot](1,2,2)
a=juwplot.wtfrshow(tfrpo2./maximum(tfrpo2,1),dt,t[1],t[end],fin[1],fin[end],0.7)
PyPlot.title("Livingston")
PyPlot.xlabel("time [s]")
Out[22]:
In [ ]: